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We present a linear stability analysis of a planar metal electrode during steady electrodeposition. 
We extend the previous work of Sundstrom and Bark by accounting for the extended space-charge 
density, which develops at the cathode once the applied voltage exceeds a few thermal voltages. In 
accordance with Chazalviel’s conjecture, the extended space-charge region is found to greatly affect 
the morphological stability of the electrode. To supplement the numerical solution of the stability 
problem, we have derived analytical expressions valid in limit of low and high voltage, respectively. 


I. INTRODUCTION 


One of the most interesting aspects of systems, involv¬ 
ing transport between matter in different phases, is their 
tendency to become morphologically unstable and de¬ 
velop ramihed growth patterns. Well known examples 
include snow flake formation and dendritic growth dur¬ 
ing metal solidification Bi- A particularly interesting 
and challenging growth problem is encountered in elec¬ 
trodeposition from an electrolyte onto an electrode 
[TTj| . Whereas snow flake formation and solidification are 
mainly driven by diffusion of water vapour and heat, re¬ 
spectively [l|, 0 , electrodeposition is driven by electroini- 
gration in addition to diffusion [T^, Ell • For this reason, 
the electrodeposition rate can be driven to exceed the 
diffusion limit, at which point the system enters a nonlin¬ 
ear regime not encountered in the purely diffusion-driven 
systems. One of the features of this nonlinear regime 
is the development of a nonequilibrium space-charge re- 
gion, which extends from the cathode into the electrolyte 
[3 El]- Already in 1990, Chazalviel realized that 
this extended space-charge region is crucial to the un¬ 
derstanding of ramified growth during electrodeposition 
E3- Nevertheless, there has been very little work which 
actually takes this effect into account. 

In this paper we investigate the morphological stability 
of the cathode during electrodeposition in both the lin¬ 
ear and the nonlinear regime. We follow the approach of 
Sundstrom and Bark |l6l |. and investigate steady elec¬ 
trodeposition in a system composed of an electrolyte 
sandwiched between two, initially planar, metal elec¬ 
trodes. We solve the stability problem numerically and 
find that the higher the applied voltage difference is, the 
more unstable the electrode surface becomes. Also, the 
most unstable wavelength becomes smaller as the voltage 
bias is increased. 

In addition to solving the stability problem numeri¬ 
cally, we derive analytical expressions for the perturba¬ 
tion growth rate, valid in the low and high voltage limit, 
respectively. In deriving these expressions, we make use 
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FIG. 1. Sketch of the studied system with lower and upper 
electrode surfaces at a; = fe{y,t) and x = 2 + /„(?/,t), re¬ 
spectively. The coordinates are given relative to the moving 
frame of reference, following the mean speed of the electrode 
surfaces, and normalized by half the electrode spacing L. 


of an accurate analytical model for the extended space- 
charge region, which we presented in a recent paper E^ ■ 


II. MODEL SYSTEM 

Following Sundstrom and Bark 0 , we consider a bi¬ 
nary electrolyte trapped between two co-planar metal 
electrodes at a; = 0 and x = 2L. The electrolyte has 
initial concentration cq and is assumed symmetric with 
valence Z. The coordinate system is moving in the neg¬ 
ative ^-direction with velocity U, which is related to the 
mean movement of the electrode. We consider the dilute 
solution limit, in which the effect of the moving coordi¬ 
nate system is negligible everywhere except in the surface 
evolution equation. A sketch of the system is shown in 
Fig.ffl 

In the analysis, we investigate the stability of y- 
dependent perturbations along the i-direction. However, 
our analysis is general and applies to perturbations along 
any direction in the yz plane. 


III. GOVERNING EQUATIONS 

The current densities of either ion are given as 
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where we have non-dimensionalized the currents J± by 
the limiting currents 2D^cq /the electrochemical po¬ 
tentials /ij_ by fcgT, the electric potential </> by the ther¬ 
mal voltage Vt = k^T/e, the coordinates by half the 
electrode spacing L, and the concentrations c± by the 
initial concentration cq. Normalizing the time by the 
diffusion time tg = L'^/{2D+), the non-dimensionalized 
ion-conservation equations become 

^dtC± = -V -J^. ( 2 ) 

At the electrodes, the current of anions vanishes, while 
the current of cations is given by a reaction expression 


Eq. (j3bl) . condition ® corresponds to ascribing the en¬ 
tire current into the upper electrode to electromigration. 

Finally, since the anions can not enter or leave the 
system the total number of anions is conserved, 

I (c_ - 1 ) dE = 0 . ( 10 ) 

JQ 

We introduce functions x = fpijj) describing the posi¬ 
tion of the upper and lower electrode u and The time 
evolution of fp is determined by the single-ion volume 
and the current into the electrode, 

{dtfe-U)ea;-ni =-a^cgne-J+, Anode, (11a) 

{dtfu -U)e^-nu = -a^cgtiu ■ J+, Cathode. ( 11 b) 


Up - J- = 0, (3a) 

'^P ' ^-t- “ ktp^ (^b) 

where Rp is the reaction rate at the lower and upper 
electrode, respectively, as indicated by the subscript p = 
£, u. We model the reaction rates Ru and Ri using the 
standard Butler-Volmer expression (l^ . 


Rp = Kg 


l^_^g-li^+aZ(<p+Vp) _ ^-'tK-{l-a)Z(<p+Vp) 


( 4 ) 


where Kg is the dimensionless version of the dimension- 
full rate constant kg for the electrode reaction, 
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kg 
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Above, Vp is the normalized electrode potential, k is 
the normalized curvature of the surface, a is the charge- 
transfer coefficient, and 7 is the non-dimensionalized ver¬ 
sion of the dimensionfull surface energy 7 , 


Here, the filling factor a?cg is much less than unity, since 
we are dealing with dilute solutions. The normalized 
velocity U of the coordinate system accounts for the mean 
current into or out off the electrodes, and dtfp accounts 
for local deviations from the mean current. 

The curvature k and the normal vectors are related to 
the surface function fp by 
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In defining the above equations and boundary condi¬ 
tions, we have chosen slightly different normalizations 
than in Ref. [l^ , the main difference being that we allow 
for a non-zero space charge density. 


IV. PERTURBATION 


0^7 

k^TL' 
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The electrostatic part of the problem is governed by 
the Poisson equation. 


2\p7‘^(j}= -p= -Zc+ + Zc-, (7) 


where the non-dimensional Debye length Ap is given as 


T _ 


with An = 


' k-gTep, 

2e^cg 


( 8 ) 


For simplicity, and to be in accordance with most previ¬ 
ous work, we choose not to explicitly model the Debye 
layers adjoining the electrodes. Instead, we apply the 
boundary conditions ([3]) just outside the Debye layer. 
Following Ref. [13 we implement the boundary condi¬ 
tion 


n„ • Vc+ = 0, (9) 

at the upper electrode, to reflect the minimum in c+ 
at the outer edge of the Debye layer. Together with 


The stability of the problem is investigated using linear 
perturbation theory. That is, we impose a small pertur¬ 
bation on a steady-state base state, and investigate how 
the perturbation evolves. The base state is identified by 
a superscript ” 0 ” and the first-order perturbation by su¬ 
perscript ” 1 ”, 


fpiy^t) - 

- fp{y,t), 
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s (j)°{x) + (j)^{x,y,t). 

(13c) 

In first-order perturbation theory, we substitute the 
second-order factor y^l -|- [dyfpY in Eq. (|T^ by unity. 

~ ^ydyfl 1 

Ttu ~ ^x ^y^y f ui 

(14a) 
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To find the field values at the perturbed surface, we Tay¬ 
lor expand to first order and obtain 

+ dx(p\ofeiy,t) + (15a) 

V<(>(/i,?/,t) edy(l)^\gey (15b) 
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Similar expressions apply for c± and at the upper elec¬ 
trode. Evaluating the reaction rate at the lower electrode 
and expanding to first order, we find 

Re ft! Re+R], 


^ - j) paZ{4P+Vi) _ -{l-a)Z{4,°+Vi) 

^ = „aZ{,l>° + Vi) 

Kn 


(16a) 
(16b) 

c+ -f d^ePfl (16c) 

c+ ( - ySy/i + olZ[( j)^ + da:(l>°fe ]) 

- 

- (1 - a)Z[(j)^ +da;(l)°fe] 

where all fields are evaluated at x = 0. Similar expres¬ 
sions apply at the upper electrode. 

Hence, the full zeroth-order problem becomes 


— e 


-{l-a)Z{^°+Vf) 


0 = -5,4, 
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24 = T Zcld^v , 
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(17a) 
(17b) 

245^4 =-44-c'i) = -p^ (17c) 

with the following boundary conditions and constraints 

4 ( 0 ) = 0 , 4(2) = 0 , (18a) 

4 ( 0 ) =-4, 4(2) = i?°, (18b) 

[ (c° - 1 ) dx = 0, 5,c+(2) = 0, (18c) 

^0 

and the mean growth velocity U derived from Eq. dD, 
U = a^coJl. (19) 

Similarly, the first-order problem is given by 


D+ 




D 




24 = -Vc^ T Zc^Vcj)^ T Zc}^V(, 


24 v 4 ^ =-Z{c\- Zc\), 


(20a) 

(20b) 

(20c) 


AD V 

and the boundary conditions, 

e,-4(2) = 0, e, •4(0) = 0, (21a) 

e,-4(2) = 4, e,-4(0) = -4, (21b) 

44(2)4 + 5,c42) = 0 , (21c) 

together with the first-order electrode growth rates dtf} 
and 5 (/m derived from Eq. (IllL 


dtfl = a^coR\, 


dtfu = -aVoi?„. 


( 22 ) 


To hnd the eigenmodes, we make the following har¬ 
monic ansatz for the first-order fields, 


c 4 x,^/,^) = ci(x)e™^ 

4(a:,y,0 = 

4(y,^) = 4eI^‘+*'=^ 


(23a) 

(23b) 

(23c) 


where E is the nondimensional growth rate of the per¬ 
turbation, and k is the wavenumber of the transverse 
eigenmode. Eor convenience we also define 


4 = Ryt+^i^y. 


(23d) 


With this ansatz, the first-order bulk equations become 


2^rci = -fc2(ci ± Zc^4) 
D± 


(24a) 


24(520*-44) = -Z(c;- 4), (24b) 

and the first-order reaction rate at the lower electrode is 

TD* 

-the 
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(25) 

yields 


Inserting the ansatz in the growth equations 

r4 = 4co4, r4 = -4coi?:. (26) 

V. ANALYTICAL RESULTS 

For large wavenumbers, fc ^ 1, we can neglect fe and 
the left hand side in Eq. (l24aE Analytical expressions 
for the growth rate can then be obtained in the limit of 
overlimiting and under limiting current, respectively. In 
Appendices lAl and iBl we find that the growth rate can be 
expressed as 


r = a'^cokJ 


i C -74 

1 + k ■ 


(27) 


With the usual assumption a = ^, simple expressions for 
1 can be obtained. In the limit of underlimiting current 
< 1 it becomes 




jo 


1-JO 


i + + ( 1 -JO) 


(28) 


We note that the underlimiting expression is nearly iden¬ 
tical to the one already derived by Sundstrom and Bark 
M- In the limit of over limiting current jO > 1, we ob¬ 
tain 




--I-’ 


(29a) 


Ar 


' 2 J 0 

1 “ 4 


(29b) 




















4 


The critical wavenumber kc, where the perturbation is 
marginally stable, is found to be 


fcc — 



(30a) 


and the wavenumber fcmax, at which the growth rate is 
maximum, is given as 


2-e7 + 2yr^^ 


2-a + 2yr^ 


-1/3 


- 1 


(30b) 


We note that the analytical model takes the zeroth- 
order current density as input variable through If 
one wants the results as a function of the potential drop 
instead, a model of the system’s current-voltage charac¬ 
teristic is needed. For simplicity, we just use the nu¬ 
merically calculated current-voltage characteristic in the 
following. 

To compute the results without reference to a numeri¬ 
cal solution, an analytical model for the system’s current- 
voltage characteristic is required. Such a model can be 
found in our previous work |15| . To obtain the total volt¬ 
age drop over the system, the interfacial voltages from 
Eq. (I16bl) should also be taken into account. 


VI. NUMERICAL SOLUTION 

The numerical simulations are carried out in the com¬ 
mercially available finite element software COMSOL 
Multiphysics ver. 4.3a. Following our previous work 
[S [H the zeroth- and first-order problems are 
rewritten in weak form and implemented in the math¬ 
ematics module of COMSOL. In the first-order prob¬ 
lem we set the parameter /„ to unity, meaning that the 
magnitude of the remaining first-order fields are given 
relative to the amplitude of the upper electrode pertur¬ 
bation. To limit the parameter space, we choose fixed. 


TABLE I. Fixed parameter values used in the numerics. 


Parameter 

Symbol 

Value 

Cation diffusivitvflSl 

D+ 

0.714 X s"^ 

Anion diffusivity[l^ 

U_ 

1.065 X 10"®m^ 

Ion valence 

Z 

2 

Surface energy 

7 

1.85 Jm"^ 

Temperature 

T 

300 K 

Permittivity of water 

€w 

6.90 X 10"“Fm-i 

Charge-transfer coefficient 

a. 

1 

2 

Reaction constant^ 

ko 

9.4 X 10^V“^ 

Diameter of a copper atom^ 

a 

0.228 nm 


^ Calculated using the exchange current /q = 30 A m~^ from 
Ref. [l^ and fco = loKZe). 

The cube root of the volume per atom in solid copper IlSl . 



FIG. 2. (Color online) Zeroth-order cation concentrations c(J. 
shown in full (black) lines and zeroth-order charge densities 
jZ shown in dashed (red) lines. The inset shows the fields 
close to the electrode. In the simulation the parameter values 
Co = 10 mM, L = 10 pm, and Vo = {1, 5,12, 30} were used. 


physically reasonable values for the parameters listed in 
Table n The values are chosen to correspond to copper 
electrodes in a copper sulfate solution. We note that the 
surface tension is quite difficult to determine experimen¬ 
tally, and most measurements are carried out at temper¬ 
atures around 1000 °C [lll,!!^. Ab initio calculations can 
give some impression of the behaviour at lower tempera¬ 
tures [ 2 ^, but these can hardly stand alone. Extrapolat¬ 
ing the linear fit of Ref. [2lj down to 0 K yields surface 
tension values close to those obtained from ab initio cal¬ 
culations in Ref. [ 2 ^. This makes it somewhat plausible 
to apply the model from Ref. [2l[ in the region of interest 
around 300 K. This yields a copper-gas surface energy of 
1.92 J/m^. The contact angle at the copper-water inter¬ 
face is very small [ 2 ^ , so finding the copper-water surface 
energy is just a matter of subtracting the surface energy 
of water from that of copper. The resulting surface en¬ 
ergy is 7 ~ 1.85 J/m^, as listed in Table H) A great deal 
of uncertainty is also associated with the value of fcoj and 
the value of a = ^ is largely a matter of convention. 

These choices leave us with three free parameters, 
which are the bias voltage Vb, the electrolyte concentra¬ 
tion Co, and the system length L. 

The solution procedure is as follows: First, the zeroth- 
order problem is solved for a given set of parameters. 
Then the first-order problem is solved for a range of 
wavenumbers k. For each k value, the corresponding 
growth rate F and perturbation amplitude of the lower 
electrode, Fe, are obtained. 

In Fig. [21 the zeroth-order cation concentrations c!|_ 
and space-charge density are shown for cq = 10 mM, 
L = 10 pm and varying bias voltage Vb- It is seen, that 
when the bias voltage exceeds Vb — 12, local electroneu¬ 
trality is violated near the cathode. For Vb = 30 the 
nonequilibrium space-charge region extends far (0.04L) 
into the electrolyte. 
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FIG. 3. (Color online) The growth rate F plotted versus the 
perturbation wavelength A for Vo = 30, co = 10 mM, and 
L = 10 iim. The full (black) line shows the growth rate 
obtained from numerical simulations, and the dashed (red) 
line shows the growth rate according to the analytical model 
Eq. (na. For perturbation wavelengths smaller than the crit¬ 
ical wavelength Ac = 51 nm the system is stable and for larger 
wavelengths it is unstable. At the most unstable wavelength 
Amax = 110 nm the growth rate is Fmax = 0.0193. 


A. Results 


For plotting purposes we introduce the dimension- 
full perturbation wavelength A = 27rL/fc. In Fig. [3l 
the growth rate F is plotted versus A for Vq = 30, 
Co = 10 mM, and L = 10 pm. Visible in the figure is 
a stable region for wavelengths smaller than the critical 
wavelength Ac = 51 nm, and an unstable region for larger 
wavelengths. The most unstable wavelength we denote 
Amax, and the corresponding growth rate we denote Fmax- 

To enable a more compact representation of the data, 
we introduce a gray-scale contour plot of the magnitude 
of F, as illustrated in Fig. IH Here, F is plotted versus 
the wavelength A for Vq = {5, 10, 15, 20, 25, 30}. The 
gray scale in the A-Vb plane is created by projecting the 
F values from the above curves onto the plane. The solid 
(blue) line in the (A, Vb)-plane marks the crest of the hill, 
thus representing the most unstable wavelength for each 
value of Vo . 

In Fig. [SJ we make use of the contour plots to show 
results for twelve sets of (cq, T)-values. In each contour 
plot, F is normalized by its maximum value, which is 
given above each plot. Shown in thick lines are Amax 
in bright (yellow) and Ac in black. The corresponding 
analytical results are shown in dashed (blue) and dotted 
(green) lines, respectively. The thin black lines show con¬ 
tours, where F equals {0.01,0.2,0.7} times the maximum 
value. There is a clear tendency in all of the panels that 
the growth rate F increases rapidly with Vq , and the most 
unstable wavelength decreases as Vq increases. Across the 
panels, the maximum growth rate is seen to increase for 
increasing cq and increasing L. Also, the most unstable 



FIG. 4. (Color online) The growth rate F plotted versus the 
perturbation wavelength A and voltage Vo for co = 10 mM, 
and L = 10 pm. The (cyan) space curves are plots of F versus 
A for Vb = {5, 10, 15, 20, 25, 30}. The shade of the in plane 
contour plot is based on the logarithm of F, which is why there 
are no contours in the low A limit where F is negative. The 
thick (blue) in plane line marks the crest of the hill, i.e. it 
marks the most unstable wavelength for each value of Vb- 


wavelength Amax and the critical wavelength Ac become 
smaller as cq increases and as L decreases. 

A common feature seen in all of the panels, is the 
kink in the Vb-versus-Amax and Vb-versus-Ac lines. At 
this kink, the slope of the lines changes markedly. The 
kink is located at the voltage, where the current reaches 
the limiting current, and it thus signifies that there is a 
qualitatively different behavior for over- and underlimit¬ 
ing current. This qualitative difference between the two 
regimes is in accordance with the analytical models. We 
also see that the kink voltage changes with cq and L. 
Specifically, it increases with cq and decreases with L. 
The main reason for this behavior is easily understood 
with reference to the zeroth-order Butler-Volmer reac¬ 
tion expression (Il6bl) . Setting the current in the system 
to the limiting current J° = 1, the reaction rates at the 
electrodes become 


- TUq 


^_^^aZ(<fi+Vp) _ ^-(l-o,)Z(4,+Vp) 


(31) 


At the cathode, the first term in the bracket dominates, 
and at the anode the other. Therefore, both potential 
drops over the electrode interfaces scale as 

AF^-ln(ifo) = ln(^^^^ 

which increases monotonically with increasing cq/T. As 
a consequence, the total potential drop at the limiting 
current also increases with increasing co/T, just as ob¬ 
served in Fig. m 

In addition to the instability growth rate F, which gives 
a time scale for the development of instabilities, it is use¬ 
ful to have a measure for the characteristic instability 
length scale. For instance, we would like to estimate the 
thickness of the deposited layer, when instabilities start 
to develop. We define this instability length scale as the 
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FIG. 5. (Color online) Contour plots of F plotted versus wavelength A and voltage Vo for co = {1 mM, 10 mM, 100 mM} and 
L — {1 pm, 10 pm, 100 pm, 1 mm}. In each plot, F is normalized by its maximum value, and the contours are logarithmically 
spaced. The maximum value Fmax of F is given on top of each plot, and the point where the maximum value is attained is 
indicated with a dark (red) circle. The three thin black lines in each plot indicate contours where F equals 0.01, 0.2, and 
0.7 times Fmax- The thick bright (yellow) line marks Amax for each value of Vb, and the dashed (blue) lines mark the two 
corresponding analytical limits. The thick black line marks Ac for each value of Vo, and the dotted (green) lines mark the two 
corresponding analytical limits. 


product of the zeroth-order growth rate Eq. (1191) and the 
instability time scale at the most unstable wavelength 


Ep = L 


a^co J4. 

^ max 


(33) 


where the pre-factor L ensures a dimensionfull expres¬ 
sion. In Fig. [6l we plot the instability length Lp ver¬ 
sus applied voltage Vb for L = 100 pm and varying cq. 
The most unstable wavelength Amax is also plotted in the 
same figure (dashed lines). It is seen that Lp decreases 
as Vb increases, but for small voltages Lp is largest for 
high concentrations, while the opposite is true for high 
voltages. The reason for this reversal is that the inter¬ 
facial voltage drops are largest for large cq. At small 
voltages the bulk driving force in the systems with large 
Co is therefore small, and this causes the system to be 
less unstable than the low cq systems. We also see that 
Amax scales in the same way as Lp. While the reason for 


this is not immediately obvious, it is seen to follow from 
the analytical expressions. Inserting Eq. (E71) in Eq. (1551) 
yields 


Lp = 




2ttL 


27T 




2 ’ 


(34) 


which confirms the approximate scaling between Lp and 
Amax. The connection between Lp and Amax implies that 
Amax sets the scale, not only for the variations in the 
horizontal direction, but also for variations in the verti¬ 
cal direction. We might therefore expect that the rami¬ 
fied electrodeposits, emerging at much longer times than 
r“ax: have a universal length scale roughly set by Amax- 
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FIG. 6. (Color online) The instability length scale Lt (full 
line) and most unstable wavelength Amax (dashed line) plotted 
versus bias voltage Vo- The concentration varies between the 
values Co = {1 mM, 10 mM, 100 mM} and the length L = 
100 pm was used. 


VII. DISCUSSION 

The main feature, which sets our work apart from pre¬ 
vious stability analyses of electrodeposition, is the inclu¬ 
sion of the overlimiting regime. Presumably, this regime 
has so far been avoided due to the non-linearities arising 
at overlimiting current, which necessitate a more com¬ 
plicated treatment. However, the overlimiting regime 
is highly relevant for ramified growth problems [7|, y. 
As seen in Fig. [SJ the instability growth rate increases 
markedly in the overlimiting regime, and there is also a 
change in qualitative behavior between the two regimes. 
Of course, the conclusions we reach, based on our model, 
are only strictly valid for planar electrodes. It does, how¬ 
ever, seem reasonable to expect that the most unstable 
wavelength Amax is comparable to the characteristic di¬ 
mensions encountered in a ramified growth experiment. 
Our analysis can thus be used to rationalize experimental 
results. In this regard, our analytical models are particu¬ 
larly useful, since they allow for easy computation of the 
key parameters for other systems than the one treated 
here. 

Perhaps the most important application of the stabil¬ 
ity analysis, is as a means of validating more elaborate 
numerical models of ramified growth. A model of rami¬ 
fied growth must necessarily deal with a moving interface 
and this, as well as other complications, make for highly 
complex numerical models. To validate such models it 
is very useful to have a comparatively simple model, like 
the present one, to benchmark against in the relevant 
limit. Indeed, this was what originally motivated us to 
treat the stability problem. 

An obvious shortcoming of the given analysis, is the 
restriction to a steady-state zeroth-order solution. The 
principal reason for this choice is that it makes for a 
simpler problem. Furthermore, the numerical ramified 
growth model, to which we wish to compare our model, 
is at present also restricted to quasi-steady state. In time, 


we wish to extend both models to the fully transient 
regime. 

There is, however, some physical justification for mak¬ 
ing the steady-state assumption. As seen in Fig. O the 
growth rate F is considerably smaller than unity in a 
large part of the investigated parameter space. The time 
it takes the system to reach steady state is given by the 
diffusive time, which in our normalization has the value 
one. Thus, as long as F is much smaller than unity, the 
system reaches steady state long before any instabilities 
build up. In this case it is therefore justified to assume 
steady state. It should be noted that in this argument 
we make the reasonable assumption that the true growth 
rate in the transient regime does not significantly exceed 
the steady-state value. 

To model the reaction rate at the electrodes we use 
the standard Butler-Volmer expression, which stands at 
the core of much electrochemistry. The Butler-Volmer 
model is, however, largely a phenomenological model, 
that does not necessarily apply in all regimes [1^. As 
applied here, the model does not distinguish between the 
potential drop over the Debye layer and the potential 
drop in the narrow interface region. This does not seem 
quite right, but it could be fixed with relative ease by im¬ 
plementing the Frumkin correction to the Butler-Volmer 
model. Nevertheless, to be in accordance with most pre¬ 
vious studies, and because the Frumkin correction would 
introduce another unknown parameter, we have chosen 
not to make this correction. 


VIII. CONCLUSION 

We have successfully solved the stability problem in the 
under- and overlimiting regime for the case of a copper 
sulfate solution trapped between two copper electrodes. 
In addition to the numerical solution of this particular 
problem, we have derived analytical solutions valid in 
either the overlimiting or the under limiting limit. The 
behavior in the overlimiting regime differs qualitatively 
from the behavior in the underlimiting regime, and we 
find that the electrode becomes increasingly unstable as 
the current is increased above the limiting current. The 
stability analysis, and in particular the analytical limits, 
are valuable both for rationalizing experimental results 
and for validating more elaborate numerical models of 
ramified growth. 

Appendix A: The electroneutral limit 

In the limit where the electrolyte is locally electroneu¬ 
tral and the time derivatives in the first-order transport 
problem are negligible, analytical solutions to the prob¬ 
lem can be obtained. Setting the point of zero electro¬ 
static potential at a; = 1, it is easily found that 

c = c+ = c_ = . (Al) 





It follows that c = +e^^°Z(j?-, and thus 


and solving for C, we obtain 


c° = and Zcj)^. (A2) 

Solving the zeroth-order problem yields 

cO = l-jO(x-l), Z(jP =\n{l- J°{x-l)). (A3) 

Using the electroneutrality assumption in Eq. (I24ap we 
find 


Q = dlc*-k^c*. (A4) 

This equation has two solutions, but as long as the per¬ 
turbation wavelength is considerably smaller than the 
electrode spacing, the solution which increases with x 
is dominant 


(^gfc(x- 2 ) 


(AS) 


where C is a constant to be determined. From Eq. (ESI) 
we then find 


1 c* C 

^ 'Zl- J^{x-iy 


(A6) 


At the upper electrode, a; = 2, the first-order reaction 
rate is (we set = 1) 


T>* 

_ 

Ko ~ 




— e 


-{l-a)Z{4>° + V^) 


-{l-a)Z[cj)* +d^cjy] . (A7) 


We can simplify this expression by using 


^-{l-a)Z(4P+V^) ^ 


K 

if o’ 


(A8) 


from the zeroth-order reaction expression. The first- 
order reaction expression then becomes 


-2Xoe“^(-^“+^“) -k + {l- a)j^ 


(A12) 


The growth rate can be expressed as 

r = -a^coJ* = a^cokC, (A13) 


SO we have 


r = o/cokT 


,2ifoe“^(^”+^“) - 7fc' - (1 - a)^ 
2Koe^zic^°+V-^) +k-{1-a) 


(A14) 


For the special case where a = ^ we find from the zeroth- 
order current that 


2Koe°‘ 


2(0“-I-14) 


JO 



(A15) 


Inserting this yields 


jO 

1-J“ 

a+-\/l + 4(^)'(l-JO) 

— 

JO 

1-JO 

a+^l + 4(^)"(l-JO) 

+ k 


(A16) 

an expression which accurately replicates the numerical 
results at low voltages and low Ap. 

To test whether the time derivatives in the first- 
order problem really are negligible, we compare the time 
derivative term 2Vdy with the transverse diffusion term 
k'^cy. Since Eq. (IAI6|) implies T < a^cokJ^, our assump¬ 
tion is justified if 


2a^CoJ° < k. (A17) 


TD* 

Kr 


= e 


Kn 


c* +d,c° + clz{(j)* + 


^k'^-{l-a)Z{cj}*+d,ct)°) 


(A9) 


Evaluating the fields at x = 2, this expression becomes 


Consequently, because a^co 1 for dilute systems and 
jo is of order unity, it is justified to neglect the time 
derivative, unless the perturbation wavelength is much 
larger than the electrode spacing. 

The critical wavenumber kc is found by setting the 
nominator in Eq. (lAlBp equal to zero. 


:^ = 2(C-jO)e“^(^“+^“) 

A'o 


+ 


ifn 


— (I — a) 


C- jo 


1-JO 


(AlO) 


The current into the upper electrode is J* = —dxC* = 
— kC, meaning that 


-kC = K = 2{C - jO)iFoe“^^‘^“+'"“^ 

-f jOyfc^ - (1 - a)Jo ^_ , (All) 


kc = 



(A18) 


where we have introduced the parameter 




JO 


1-JO 


« + i/i+d$l d-J") 


(A19) 


To find the wavenumber /cmaxi at which T attains its max¬ 
imum Tmax, we set the derivative of T equal to zero and 
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solve for fc, 


region (ESC) they take the simple forms 


k — 

^max — 




1/3 




(A20) 


with the asymptotic solutions, 


Zj 

72 Ad 


-I -1/2 


X —1 — 


JO 




X —1 — 


dx(j>°{x) R! -^Tjo 
Ad 


1/2 


-3/2 


X — 1 — 


JO 


(B6) 

(B7) 

(B8) 


I (^) ^ , for 7^ » 1, 

I (ft) ' - i foi 7^ < 1- 


(A21) 


Appendix B: The strongly non-linear limit 


In the limit where the driving force is very large, some 
of the terms in Eqs. (I24al) and (I24bl) become dominant, 
which makes an analytical solution of the problem possi¬ 
ble. 

If the system is strongly driven, the field gradients are 
large close to the upper electrode, and this makes the 
electrode surface much more unstable. It follows that a 
larger k value is needed for the surface tension to stabilize 
the system, so the most unstable value of k will be larger 
than for less driven systems. In the strongly driven limit, 
we might therefore expect that Eq. (I24bl) largely is a 
balance between d'^cj)* and k'^cj)* in the region of interest. 
This leads us to making the ansatz 

(j)* = (Bl) 

where $ is a constant. We now consider Eq. (I24al) for 
the cation concentration, neglecting the left hand side 


The width of the ESC is given as Tesc = I — 1/J°, so in 
the region close to the electrode, compared to the width 
of the ESC, the fields can be written as 



Evaluating we find 


OxC 


+ 


Z 3c° 


6 * - Z- 


2 Tesc 


(B9) 

(BIO) 

(Bll) 


(B12) 


which is seen to be much smaller than Zc^dx4>* if 


2fc>--^, (B13) 

TeSC 

that is, if the perturbation wavelength satisfies 


A <C 47rLESC- (B14) 


Similarly, we find that Zc*^dx4>^ is much smaller than 
Zc%dx(t>* if 


0 = - dxc\ - Zcldx4P - ZclOx^^ 

-7(c;-bZc°. (/)*). (B2) 

We assume that the terms dxc'^ and Zc\dx(j)^ are negligi¬ 
ble compared to Zci^dx(t>* and insert the ansatz Eq. (IBII) 

0 R Zdxc%k(l)* + Zc%k‘^(j)* - k^{cX + Zc%(t)*) (B3) 

K. Zdxc%k(t)* -k^cX, (B4) 

implying that 

c; R jdxclc^*. (B5) 

To test the assumptions leading to this result, we need 
expressions for c(|_, dxc\ and dx4>^■ Erom Ref. [T^ we 
have such expressions, and in the extended space-charge 


A^ < 


8tt^ Ad 

^ cl{2)' 


(B15) 


Finally, the ansatz Eq. (IBII) is justihed if 2\^k‘^(f>* ^ 
Zc*^_^ which is equivalent to 


This last requirement is seen to follow if the two first 
requirements Eqs. (IBI4I1 and (IBI5I1 are fulfilled. 

In the strongly driven regime, where Eqs. (IB14I) 
and (jBlhp are satisfied, the first-order current is approx¬ 
imately 


2j; R -Zc%dx(t>* = $, (B17) 

at the upper electrode. The zeroth-order diffusive contri¬ 
bution is also very small at the upper electrode, meaning 
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that we can simplify Eq. (IA9I) 


i?: 




K, 


- (1 - a)Z{(/)* + d^(l)°) 


Koe' 


+ K 


,aZ(0O+V„ 


- 2r 


7 r-(l-a) 


2J°\ 

'-tJj 


Inserting i?* = Ri — iZfcc° $ we find 


(B18) 


(B19) 


2^ggaZ(0'>+V7) _ _ (1 _ 


_ju„o * _ )u 70_ 

2 + 2Xoe“^('^“+^“)+fc-(l-a)^ 


(B20) 


and since E = — o^cq J* 


- 7P _ (1 _ 

r = a3cofcJ°-77^-:^77;r^7^^- 7 -(B21) 


2Koe°‘^(^°+^-') +k-{l-a) 


= a^cokJ^ 


,7“ 

2a-l+sJ\ + A{^)\\ 

— 

JO 

IF 

0 + 

2a-1 + ^l + 4(^)%^ 

+ k 


(B22) 


where the second expression assumes a = i, so that 
Eq. (jA15l) is valid. Like in the electroneutral limit ne¬ 
glecting the time derivative in the first-order problem 
is justified, unless the perturbation wavelength is much 
larger than the electrode spacing. The expressions (IA18I) 
and (IA20|) are also valid for the strongly nonlinear limit, 
if instead of Eq. (IA19|) we use 




jo 


2cr — It 



(B23) 
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